Post-pubertal developmental trajectories of laryngeal shape and size in humans

Laryngeal morphotypes have been hypothesized related to both phonation and to laryngeal pathologies. Morphotypes have not been validated or demonstrated quantitatively and sources of shape and size variation are incompletely understood but are critical for the explanation of behavioral changes (e.g., changes of physical properties of a voice) and for therapeutic approaches to the larynx. This is the first study to take this crucial step and results are likely to have implications for surgeons and speech language pathologists. A stratified human sample was interrogated for phenotypic variation of the vocal organ. First, computed tomography image stacks were used to generate three-dimensional reconstructions of the thyroid cartilage. Then cartilage shapes were quantified using multivariate statistical analysis of high dimensional shape data from margins and surfaces of the thyroid cartilage. The effects of sex, age, body mass index (BMI) and body height on size and shape differences were analyzed. We found that sex, age, BMI and the age–sex interaction showed significant effects on the mixed sex sample. Among males, only age showed a strong effect. The thyroid cartilage increased in overall size, and the angulation between left and right lamina decreased in older males. Age, BMI and the age–height interaction were statistically significant factors within females. The angulation between left and right lamina increased in older females and was smaller in females with greater BMI. A cluster analysis confirmed the strong age effect on larynx shape in males and a complex interaction between the age, BMI and height variables in the female sample. The investigation demonstrated that age and BMI, two risk factors in a range of clinical conditions, are associated with shape and size variation of the human larynx. The effects influence shape differently in female and male larynges. The male–female shape dichotomy is partly size-dependent but predominantly size-independent.


Results
Shape variation. The first two shape axes of the mixed sex sample accounted for 31.7% and 10.9% of the explained variation, respectively (Fig. 1A) and regressions of principal components (PC) PC 1 and 2 scores on sex are significant (Table 1). Variation among individuals along the first principal component (PC 1) mostly reflects male-female differences (Fig. 1A). While the first PC reflects size differences (regression of PC 1 scores on centroid size: R 2 = 0.40, p < 0.001), this is driven by the fact that males are larger than females on average; the regressions of PC scores on cartilage size within each of the sexes are statistically significant for females (R 2 = 0.047, p = 0.026) but not for males (R 2 = 0.009, p = 0.33) and do not explain much of the variation even in females. In other words, individuals with a larger thyroid cartilage do not differ from individuals of the same sex with a smaller cartilage in the same direction that males differ from females. The primary axis of shape variation in the thyroid cartilage relates to the angulation of the right and left thyroid laminae, the orientation of the superior border of the thyroid cartilage and the relative height of both laminae (Fig. 1B). Lower scoring individuals, which are predominantly female, have a thyroid cartilage that is anterior-posterior (AP) compressed and relatively wide superiorly with a superior border that is angled upward ventrally (Fig. 1B). The second PC of thyroid cartilage shape contrasts individuals with an AP compressed cartilage and a long superior horn and those with a more AP expanded cartilage and a shorter superior horn (Fig. 1C). Unlike variation along PC 1, variation along PC 2 does not correspond to male-female differences (Fig. 1A).
The correlation among sex, age, BMI, body height and cartilage size. Three independent variables, age, sex and BMI, were uncorrelated with each other ( Fig. 2A-C). The cartilage was larger on average in males than females (Fig. 2D). There was no correlation between age and cartilage size (Fig. 2E) and no correlation between BMI and cartilage size when both sexes were combined (Fig. 2F). Unexpectedly, centroid size was positively correlated with age in males (Fig. 2E). Average centroid size was about 4% larger in the 78 to 88-year age group than in the 18 to 28-year age group, suggesting a growth rate of about 0.67% per decade. Females showed no correlation between centroid size and age (Fig. 2E).
Males were taller than females (Fig. 2G). While body height was weakly associated with age and BMI (Fig. 2H,I), body height, as expected, showed a strong correlation with cartilage size in the mixed sex sample; the relationship between height and cartilage size is stronger in males than females (Fig. 2J). Body height explained neither the increase in larynx size in aging men (partial correlation controlling for height r = 0.428, p < 0.001), The relationship of larynx shape to sex, age, BMI and centroid size. The full MANOVA model including age, sex and BMI as main effects, and their two-way interactions, accounted for 20.3% of shape variation in thyroid cartilage shape (19.2% when the non-significant interactions were dropped from the model) ( Table 2). There was a significant interaction between sex and age, which we explore in greater detail below, but it accounted for only 0.8% of the shape variation. The most significant factor affecting thyroid shape was sex, which accounted for 14.9% of variation independent of the influence of age and BMI.
Sex specific variation in cartilage shape. Next, we evaluated the shape variation in the thyroid cartilage for each sex. The primary axis relates to the angulation of the right and left thyroid laminae and the orientation of the superior border of the thyroid cartilage within both males and females. The angulation was larger and the superior border points cranially in the lower scoring individuals in both sexes (Fig. 3A,B). High scoring individuals in each sex have a thyroid cartilage that is more AP compressed and relatively wider superiorly with a superior border that is angled upward ventrally (Fig. 3A,B). Lower scoring individuals in each sex demonstrate tighter angulation between the left and right thyroid lamina. In males, the lamina height was greater in high scoring individuals. The second PC of thyroid cartilage shape captures variation in the relative length of the superior horn in both sexes (Fig. 3A,B). In females, the second axis contrasts also differences in relative AP distance (Fig. 3A). Next, we evaluated the effect of age, BMI, centroid size and height on cartilage shape within each sex separately. We used the same MANOVA model described above but removed the sex variable and added the height variable, with pairwise interactions again included in the model. The full model accounted for 11.2% and 8.4% of shape variation in the female-and male-only analyses, respectively ( Table 3). The interaction between age and height, and the main effects of age and BMI were significant in females accounting for 1.9%, 4.2% and 1.9%, respectively (Table 3). In males, only age was a significant factor, explaining 3.3% of shape differences.
We illustrated shape variation associated with increasing age for each sex in Fig. 3C-E since the interaction between sex and age was significant in the mixed-sex model, and age played a role in both sex-specific models. In older females, the angulation between left and right thyroid lamina increases (i.e. thyroid cartilage is AP Figure 1. Principal component analysis of shape coordinates (PC 1 and PC 2) from landmarks. Each subject is represented by a point in the two-dimensional space of PC 1 and PC 2 (A). Warped cartilage images showing shape differences in the thyroid cartilage for the first (B) and second shape axis (C). Images represent shape differences associated with the minimum and maximum extents of PC 1 and PC 2. Note a considerable overlap between the sexes on the first shape axis and a complete overlap on the second. www.nature.com/scientificreports/ compressed) and the superior border becomes angled upward ventrally (Fig. 3C). The thyroid cartilage is also more AP compressed in females with a smaller BMI, i.e. it appears as if older age and larger BMI have opposite effects on the angulation in females (Fig. 3D). We found that the angulation between left and right thyroid lamina decreases in older males (Fig. 3E) and the laminae are relatively more expanded cranio-caudally (Fig. 3E). The effect of BMI appears to be similar in males and females, however note that BMI only reached statistical significance in females (Table 3). In order to translate statistical findings into functionally relevant information, we investigated seven measures (five linear measures and two angles) that were previously identified as biomechanically relevant (Fig. 4A). Those seven measures were taken from the warped extremes for age, BMI, centroid size and height. Note that only age was significant for males and only age and BMI were significant for females (Table 3). Figure 4B illustrates direction and magnitude of the two effects on each measure. Cluster analysis. Next, we explored different morphotypes first within the mixed-sex sample and then among females and males. Using k-means cluster analysis, the 2-cluster solution for the mixed-sex sample separated most males (75.3%) into Cluster 1 and most females (87.1%) into Cluster 2 (k-means clustering; silhouette coefficient = 0.23; 20 shape axes) (Fig. 5A). Individuals in Cluster 1 were taller than in Cluster 2, and their larynx  www.nature.com/scientificreports/ size was larger (consistent with typical sex differences) ( Table 4). There were no differences in age and BMI between the two clusters ( Table 4). Next we investigated characteristics of the males and females in both clusters. Males that clustered with females in Cluster 2 (26 out of 105) were significantly younger, of shorter stature and had smaller thyroid cartilages than that of those males in Cluster 1. Indeed, Cluster 1 could be more accurately described as containing the larger individuals of each sex and Cluster 2 includes smaller individuals of each sex, even though the size difference did not reach significance within females ( Table 5). The difference between the average thyroid cartilage shapes in the two clusters ( Fig. 5D) resembles the previously described sex differences (Fig. 1C). In Cluster 2, individuals have a thyroid cartilage that is AP compressed and relatively wide superiorly with a superior border that is angled upward ventrally (Fig. 5D).
When the k-clustering analysis was performed on females only, there were four female clusters recovered that differed in age, BMI, and thyroid size (Fig. 5B, Table 6). Females in Cluster 4 were 46.4 years old on average, showed the highest BMI but smallest centroid size. Cluster 3 is the smallest cluster with only 9 individuals. This cluster had the highest average age (65.3 years) of the 4 clusters and a low average BMI (together with Cluster 1). Females in Clusters 1 and 2 were intermediate in age (59.5 and 55.6 years on average, respectively) and females in Cluster 1 had a low average BMI. Figure 5E illustrates shape differences between the 4 female clusters. The for the first and second shape axis of the single sex PCAs. In both females and males, individuals on the first shape axis are differentiated by the angulation between the left and right lamina. Images in C through J illustrate thin-plate-spline warps of regression models for age (C,E) and BMI (D,F) on a mean female and male surface rendering, respectively. The procedure first selects the individual that is closest to the group (e.g., female, male) average. Then the group average was warped toward the two ends of the slope (regression line) defined by the regression coefficients of the two models. Shape differences were exaggerated by a factor of 2 to make them more visually prominent. Note that age and BMI were significant factors for females but only the model for age was significant for males (    27 . Five linear measures and two angles were measured in the current data set using regression models for age and BMI (illustrated in Fig. 3). (B) Arrow direction indicates the direction of change of the measure with an increase in age and BMI, respectively. Arrow thickness indicates the magnitude of change between the 18-to-28-years age group and the 78-to-88years age group. www.nature.com/scientificreports/ angulation between left and right lamina is wider, and the degree of AP compression is more prominent in Cluster 3, which contains individuals with the highest average age. Three male clusters differed significantly in age and centroid size (Fig. 5C, Table 7). The youngest individuals with the smallest larynx were in Cluster 1. Individuals in Cluster 3 were oldest and their centroid size was largest. Age and centroid size of individuals in Cluster 2 were intermediate between Clusters 1 and 3. The angulation between left and right lamina was highest in Cluster 1 and lowest in Cluster 3. Interestingly the oldest/largest and smallest/youngest are not the most different in shape. Figure 5F illustrates shape differences between the 3 male clusters. The most different are the young/small vs. intermediate individuals. The angulation between left and right lamina is wider in clusters 1 and 3. The height of the laminae is taller in Clusters 2 and 3.

Discussion
A strength of this study is the use of a large stratified sample that allows for the statistical separation of multiple variables. The analysis of thyroid cartilage confirmed sex as a major source of size and shape variation. Organ size was significantly different between sexes but did not explain shape variation among the sexes, instead our  www.nature.com/scientificreports/ results support the hypothesis that additional sources of shape variation, for example age and BMI, exist within sexes. Four findings especially inform the understanding of postpubertal shape changes of the human larynx. First, the angulation between left and right lamina is a major shape difference within each sex, but also separates the average male and female shape. This does not appear to be a simple issue of size-shape (i.e., allometric) scaling. Second, the male larynx demonstrates a life-long size increase, on the order of ~ 0.7% per decade. Third, age has a significant effect on thyroid cartilage shape within males and females, as do BMI and the age-height interaction in females only. Post-pubertal changes in larynx shape differ between males and females and the shape changes associated with sex and aging have biomechanical consequences. Fourth, cluster analyses identified at least four morphotypes among females that correspond to differences in age, BMI and centroid size. In males, cluster analyses identified three morphotypes which differ by age and larynx size. Next, we discuss the relevance of the four findings and hypothesize that life-long post-pubertal trajectories of laryngeal shape changes are causally affected by BMI and body size and are, at least in part, contributing to functional changes of the larynx.
Patterns in the mixed-sex sample. The variables examined in the current study (sex, age, BMI and their respective interactions) explain about a fifth of the total shape variation in the thyroid cartilage. Even taking BMI and age into account, 15% of shape variation in the mixed-sex sample is attributable to sex. The BMI, age and height variables account for about 10% of total variation (11.2% and 8.4%, respectively in the female and male samples) in the sex-specific MANOVA models. Other genetic and environmental factors explain the remaining 80% of the shape variation.
The shape differences along the first principal component axis were associated with variation in sex. The shape difference concerns the well-known angulation between left and right lamina and the height of the laminae 12 . Table 5. Comparison of mean values (standard deviation) for age, BMI, thyroid CS and body height for males and females classified into Clusters 1 and 2 based on the mixed sex sample. The majority of subjects in Cluster 2 were male and the majority in Cluster 1 were female. The average (standard deviation) of each metric in Cluster 1 and 2 was compared between the incorrectly clustered group for males and females separately. Twosample t-tests suggest that males which were classified into cluster 2, tended to be significantly younger, smaller and with smaller centroid size. For females, none of the variables showed a difference between Cluster 1 and 2 averages. Sigific. codes: < 0.05 '*' .  Table 7. Mean values (standard deviation) for the age, BMI, thyroid size and body height variables for males in different clusters. P values indicate that the means between clusters were significantly different for age and centroid size but not for BMI and body height. Signif. codes: < 0.001 '***' . www.nature.com/scientificreports/ Males have relatively longer (antero-posteriorly), taller and more acutely angled thyroid cartilages (Fig. 3A,B). Interestingly, the primary axis of shape variation within the male and female samples shares similar shape differences to that seen in the mixed-sex sample. On the surface this appears to be a case of size-shape scaling (i.e., allometric shape differences) given that all males are, with few exceptions, larger than any female in our sample. However, this scaling effect does not hold within either sex: individuals with larger cartilages do not score higher on PC 1 in the mixed-sex analysis than those with smaller cartilages within either sex. This hints at factors other than simple scaling driving the sex-differences.

Cluster 1 (n = 37) Cluster 2 (n = 37) Cluster 3 (n = 31) P-value
Sex-specific morphotypes. Separate cluster analyses among females and among males identified different patterns of clustering. The approach appears to illustrate different laryngeal morphotypes, i.e. groups of individuals share large-scale shape features. The existence of laryngeal morphotypes among humans has been hypothesized in the context of structural adaptations for specific voice type in singing. Each voice type (e.g., soprano, alto, tenor, bass) corresponds to a pitch range which is associated with characteristics in laryngeal structure 5,6,47-52 . Among the non-professional voice users, phenotypic variation of the larynx structure has remained unexplained. The current study informs also the debate of the role of vocal organ size in determining voice features. Voice fundamental frequency is sometimes simplistically described to be a function of larynx size but should take vocal fold characteristics into account 37 . The size of the cartilaginous framework can be misleading. Thyroid cartilage shape features (e.g., relative antero-posterior extension and angulation between the thyroid laminae) that affect the shape of the vocal folds, and therefore fundamental frequency, are size-independent within sexes. Although both sexes show size differences in the larynx related to selection for larger body size, selection on males for larger larynx size related to the production of lower fundamental frequency voices 53 does not result in similar changes in the female larynx. In other words, the difference we see in the shape between sexes has evolved in the two sexes through sexual selection rather than representing an extension of a broader size-dependent scaling of shape. Both the larger size and relatively longer antero-posterior shape of the cartilage contribute to longer vocal folds documented in males compared to female 13 which contributes to the typically deeper voice in males. However, vocal fold shape, thickness and composition contribute as well. Moreover, the taller and narrower thyroid cartilage in males contributes to a longer and narrower epilaryngeal tube. The epilaryngeal tube plays an important role for vocal tract filtering 54 .
Shape changes of the thyroid cartilage could contribute to age-related changes in swallowing and breathing function of the larynx [55][56][57] . Males (relative to females), women with higher BMI and older males all have relatively more acute angles between the left and right thyroid laminae thereby narrowing the airway (Fig. 4). Shape of the cartilaginous framework is associated with remodeling of the soft tissue inside and above the larynx in a mouse model 58 . The remodeling of the soft tissue in and above the larynx could have consequences for airway patency and glottal sufficiency.
The lifelong growth of the male larynx. We found that the overall size of the male thyroid cartilage increased at a rate of about 0.7% per decade. Centroid size of the female thyroid cartilage did not increase in a systematic age-dependent manner. As far as we know, a post-pubertal size increase of the larynx has not been reported in previous studies. In general, all soft and hard tissues of the human head and neck region continue to change in shape or size after puberty 59,60 . Specifically, two other cranial cartilaginous structures increase in size throughout life similar to the male larynx. The human outer ear demonstrates sexual size differences and a lifelong size increase in both sexes [61][62][63] . Likewise, the external nose is also sexually dimorphic and its dimensions increase in both sexes up to an age of at least 65 years 64 . Testosterone sensitivity is the main driver of the laryngeal growth spurt during puberty. It remains to be seen whether an individual's hormone profile continues to affect larynx shape 65 or whether other mechanisms are responsible for the post-puberty continued growth and shape changes.
Post-puberty trajectories of shape change. Sexual dimorphism of the larynx develops during adolescence between 12 and 18 years of age 13 and maturational growth of the human body ends between 18 and 25 years of age 66 . Results from the current study suggest that laryngeal shape and size continue to change beyond the adolescent period. Those changes are different in females and males, and various environmental factors might be responsible. The post-pubertal trajectory of the male larynx reflects shape changes related to scaling. In contrast, post-pubertal shape changes of the female larynx seem more complex. Multiple factors play a role, with several potential trajectories.
In general, shape and size changes of the human body may manifest through spurts of remodeling 67 . The laryngeal shape change in adolescent males during puberty, for example, occurs relatively fast 13 and can be associated with notable voice problems (i.e., breaking voice) for extended periods of times before the adult normal male voice has manifested 68 .
A second period of hormonal changes occurs between 45 and 55 years of age associated with permanent cessation of ovarian function (menopause) in females and a reduction of testis function in males 69,70 . We observed that in females and males prominent shape differences were present between one cluster with an average age of 44 and 46 years, respectively and other clusters with average ages of 55 and older (Tables 6 and 7; Fig. 5). This suggests that the age range between 45 and 55 is an important period during which quantifiable shape changes develop.

Conclusions
We identified multiple sources of laryngeal shape variation including sex, age and BMI. The sources act differently between the sexes and generate morphotypes, i.e. clusters of similar shape, within females and within males.
Interestingly, the male larynx demonstrated a lifelong size increase. Growth trajectories play a central role in life course epidemiology. They can serve as markers of prenatal or childhood development and/or for determinants of adult health outcomes. Interestingly, both, age and BMI, are linked to numerous other anatomical variation and are identified risk factors for laryngeal pathologies 71 . Two shortcomings of the present study are the missing functional data (voice, swallowing, breathing) and the absence of a mechanistic explanation for what drives shape and size changes. Mouse models, however, show interesting potential to investigate the remodeling processes. Mouse laryngeal cartilage demonstrated considerable protein turnover rates comparable to those in muscle and mucosal tissue 72 and the mouse larynx shows also lifelong shape changes 40,58 .
The current approach of 3D reconstruction of individual laryngeal elements, shape analysis through geometric morphometrics and thin-plate spline warping to extract biomechanically relevant shape variables overcomes two important challenges of laryngeal shape analysis. First, the approach is no longer dependent on cadaver measurements but can include data from very large samples of individuals and thereby control for important factors. Second, the approach is independent from anatomical landmarks and is able to deliver measurements that are relevant for biomechanical modeling purposes 27 . This study takes an important first step of integrating 3D imaging and explores functionally important variables (Fig. 4). Next, anatomical differences can be visualized and then studied through computational modeling 37,38 . The current approach promises practicability and potential data for data mining through computer-automated analysis and integration with finite element models to link shape variation with function since the analysis is relatively fast. Any new subject can be assigned to a morphotype in less than one hour after the CT scans have been taken.

Methods
Human subject sample. All methods were performed in accordance with the guidelines and regulations by Mayo Clinic and Midwestern University specified in a collaboration agreement. This retrospective study was approved by the Mayo Clinic Institutional Review Board, and the need for informed consent was waived. The radiology database was queried to identify adult patients (age ≥ 18 years) who had undergone consecutive neck CT scans for clinically-indicated purposes. The electronic medical records were reviewed for these potential study subjects, and individuals with a history of any of the following were specifically excluded (1) head and neck cancer, (2) irradiation of the head and neck region, (3) surgery involving the skull base, neck, or upper aerodigestive tract (with the exception of tonsillectomy), (4) laryngeal dysfunction documented on physical exam or symptoms suggesting a possible laryngeal disorder (e.g. hoarseness, aspiration), (5) treatment for a systemic malignancy within the past 12 months (due to potential confounding effect on weight).
To ensure broad sampling for the variables of interest, subjects were consecutively recruited in a stratified fashion based upon sex (male, female), age (18-28; 28-38; 38-48; 48-58; 58-68; 68-78; 78-88 years), and body mass index (BMI) (< 25; 25-29.9; ≥ 30 kg/m 2 ). The BMI, weight/height 2 , is correlated with body fat used by medical professionals to screen for under-or over weight status 73 . To be eligible for inclusion, subjects had to have height and weight documented within 30 days of the neck CT acquisition. Five subjects were enrolled into each of these 42 category-combinations for a total of 210 subjects.
Body height is also an important determinant of larynx size 74 , which is itself an important determinant of laryngeal function 37 . Body height was not included in the stratification criteria for subject sampling because the projected sample size would have been very large. However, the post hoc examination of the sample revealed considerable variation in body height, which allowed the investigation of its effect on larynx shape and size. We included height as variable in the investigation of variation within females and males, but excluded height from the analysis of the mixed-sex sample given its strong covariation with sex (Pearson correlation, r = 0.71, P < 0.001).
Imaging, segmentation, 3D surface renderings, landmark coordinates and centroid size. Neck CT scans were reviewed by a board-certified neuroradiologist to ensure good image quality and technical adequacy. Although performed using different multidetector helical CT scanner models, the neck CT's were all acquired in supine position using 120 kVp with an mAs that was adjusted automatically by the scanner to maintain adequate penetration in the setting of variable body habitus. As a minimum, reconstructed z-axis slice thickness had to be < 0.7 mm with an axial field of view ≤ 28 cm and a matrix size of 512. Images were rendered with a standard soft tissue kernel. Subjects were eligible for inclusion regardless of whether or not intravenous contrast was administered for the CT. All CT datasets were anonymized prior to export for shape analysis.
Reconstructed image stacks were imported into AVIZO software (v. Lite 9.0.1) and the thyroid cartilage was segmented manually. Segmentation refers to the process of generating surface renderings of anatomical structures. Derived 3D surfaces (stereolithography format) of all specimens are available on Morphobank (http:// www. morph obank. com), Project 4345.
The 3D surfaces of the thyroid cartilage were then imported to the geomorph package in the statistical software program R 3.4.4 (R Core Team 2015) 75 to perform surface landmark placement. Landmarks were collected from surface renderings in the form of three-dimensional coordinates 76 . The landmark and semilandmark protocol captured variation in homologous structures, including structural margins and surfaces such as the thyroid lamina. Our previous work showed that 100 surface semilandmarks are sufficient to capture shape variation in the thyroid cartilage to achieve comprehensive coverage 76  www.nature.com/scientificreports/ Centroid size was used as a measure of cartilage size. Centroid size is computed as the square root of the sum of squared distances of each landmark from the centroid of the cartilage landmarks, whose location is obtained by averaging the x, y and z coordinates of all landmarks 77 . Analysis. The first step in the geometric morphometric shape analysis was superimposition of the landmark configurations to remove differences in raw scale, orientation and position in the global coordinate system. This was accomplished via generalized Procrustes analysis (GPA), which extracts shape information from the raw coordinate data by translating landmark configurations to a common location, scaling configurations to unit centroid size, and rigidly rotating configurations to minimize squared distances among corresponding landmarks 78 . As the resulting shape variables exist in a non-Euclidean shape space, we projected them to a linear tangent space prior to statistical analysis. This produces a set of transformed coordinates that reflect shape differences among cartilages independent of size.
A series of principal components analyses (PCA) were then conducted on the full sample, males only and females only to reduce the dimensionality of the data and summarize the main patterns of shape variation. We evaluated only those PCs that cumulatively explained 85% of the variation. Each PC excluded from the analysis explained < 0.5% of variation. As a result, we used 22 (men and women combined), 24 (men) and 18 (women) PCs, respectively, to explore patterns of similarity.
We began by exploring the relationship among shape variables, age, BMI, centroid size and height using Pearson correlation analyses. We then used multiple analysis of variance (MANOVA) with superimposed shape coordinates as dependent variables and sex, age, BMI, and height (in sex specific analyses only) as main effects as well as their two-way interactions using the ProcD.lm function in the geomorph package for R. Type II sums-ofsquares was applied to evaluate each main effect independent of the other main effects. Statistical significance is based on randomization of residuals via permutation, an appropriate approach for high-dimensional data such as landmarks and semi landmarks 75 .
K-means clustering (R Core Team 2021) was performed on the limited subset of principal components described above to assess whether cartilages formed distinct clusters in the shape space. The clustering method is designed to create homogeneous groups of individuals (e.g., clusters) that are maximally different from other clusters. The silhouette coefficient was used to estimate the number of relevant clusters 79 using the cluster package in R 80 . Its calculation provides a measure of how similar an object is to its own cluster (cohesion) compared to other clusters (separation). The Silhouette Score was calculated as: (b − a)/max(a,b), where a is the average intra-cluster distance (average distance between each point within a cluster) and b is the average inter-cluster distance (the average distance between all cluster centroids). The Silhouette Score can range between − 1 and 1 where a high value indicates the object is well matched to its own cluster and poorly matched to neighboring clusters. Values near zero indicate that the distance between two clusters is not significant. Based on silhouette scores, the mix-sex sample was analyzed with two clusters; three clusters were indicated for the male sample and four clusters for females (Supplemental Fig. S1). We compared the average values of the age, BMI or body height variables within each cluster using t-tests or ANOVA to assess whether the observed clustering of shape axes corresponds to variation in these factors, as appropriate.

Visualizing shape change.
To illustrate what anatomical information was encoded by individual shape axes, by regression models or by individual clusters in the cluster analyses, a thin-plate spline method was used to warp the thyroid cartilage into an estimated mean shape for a set of aligned specimens. This was achieved with geomorph's warpRefMesh function 75 . The function uses the thin-plate spline method to warp a surface model of the mean shape for a set of aligned specimens, into the shape defined by a second set of landmark coordinates. For example, the latter coordinates may represent the mean of a subsample, such as males, or may be the scaled coefficients from the regression of shape on age added to the mean configuration to represent the effect of aging.
A challenge is the translation of statistical shape differences into mechanically relevant information. We extracted two angular and five linear shape measurements with biomechanical relevance for vocal organ function from warped surface reconstructions to more directly investigate the biomechanical consequence of shape differences documented in the current study.

Data availability
Anonymized DICOM data for the neck CT scans cannot be made publicly available because of the potential to create a three-dimensional rendering of the face, which is considered to be sensitive patient health information that must be protected. Derived 3D surface reconstructions (stereolithography format) of all thyroid cartilage specimens are available on Morphobank, Project 4345.